package fem2.element.crack_sda;

import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.enu.StressType;
import fem2.material.cohesive_law.CohesiveLaw;
import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.LinAlgUtilities;

public abstract class ConstantJumpSDAElement extends SDAElement {

	protected IVector alpha; // displacement jump in global coordinates

	protected double crackMeasure; // in 2d, crack measure is crack's length. in
									// 3d, crack measure is crack's area
	protected double elementMeasure;// element measure is area in 2d and volume
									// in 3d
	protected double[] normal;
	protected double[] tangent1;
	protected double[] tangent2;

	protected StressType stressType = StressType.ALGORITHMIC_STRESS;

	public ConstantJumpSDAElement(MeshPart mp, MaterialModel mm, CohesiveLaw cohesiveLaw) {
		super(mp, mm, cohesiveLaw);
	}

	public ConstantJumpSDAElement(MeshPart mp, MaterialModel mm, CohesiveLaw cohesiveLaw,
			StressType stressType) {
		this(mp, mm, cohesiveLaw);
		this.stressType = stressType;
	}

	@Override
	public IVector getCrackOpeningParameter() {
		return alpha;
	}

	/**
	 * provide the tranformation matrix by which multiply with alpha will give
	 * local jump on the discontinuity
	 * 
	 * @param s
	 * @return
	 */
	protected IMatrix computeLocalProjectionOnCrack(double... s) {
		int dim = getMeshPart().getDimension();
		IMatrix J = new Array2DMatrix(dim, dim);
		if (dim == 2) {
			J.set(0, 0, normal[0]);
			J.set(1, 0, normal[1]);
			J.set(0, 1, tangent1[0]);
			J.set(1, 1, tangent1[1]);
		} else if (dim == 3) {
			J.set(0, 0, normal[0]);
			J.set(1, 0, normal[1]);
			J.set(2, 0, normal[2]);

			J.set(0, 1, tangent1[0]);
			J.set(1, 1, tangent1[1]);
			J.set(2, 1, tangent1[2]);

			J.set(0, 2, tangent2[0]);
			J.set(1, 2, tangent2[1]);
			J.set(2, 2, tangent2[2]);
		}
		return J;
	}

	@Override
	protected double computeDeterminantOfCrackJacobian(double... s) {
		int dim = getMeshPart().getDimension();
		if (dim == 2) {
			return 0.5 * crackMeasure;
		} else if (dim == 3) {
			return 2 * crackMeasure;// since the weight for 1 Gp on 3-node
									// triangle face is 0.5
		} else {
			throw new Error("invalid dimension");
		}
	}

	@Override
	protected IVector computeCrackTraction(double... s) {
		int dim = getMeshPart().getDimension();
		IMatrix J = computeLocalProjectionOnCrack(s);
		IVector jump = new ArrayVector(dim);
		BLAM.multiply(1.0, BLAM.TRANSPOSE, J, alpha, 0.0, jump);
		IVector t = getCohesiveLaw().computeJumpTraction(this, jump);
		IVector T = new ArrayVector(dim);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, J, t, 0.0, T);
		return T;
	}

	@Override
	protected IMatrix computeCrackStiffness(double... s) {
		int dim = getMeshPart().getDimension();
		IMatrix J = computeLocalProjectionOnCrack(s);
		IVector jump = new ArrayVector(dim);
		BLAM.multiply(1.0, BLAM.TRANSPOSE, J, alpha, 0.0, jump);
		IMatrix Cjump = getCohesiveLaw().computeJumpStiffness(this, jump);
		IMatrix C = new Array2DMatrix(dim, dim);
		LinAlgUtilities.bCBt(J, Cjump, 1.0, C);
		return C;
	}

	@Override
	public void Initialize() {
		super.Initialize();
		int dim = getMeshPart().getDimension();
		elementMeasure = getMeshPart().computeVolume();
		normal = new double[dim];
		tangent1 = new double[dim];
		if (dim == 3) {
			tangent2 = new double[dim];
		}
		crackBoundary = new ArrayList<double[]>();
	}
}
